#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

/* compile with option -lm */

/* gaussian: same alogrithm as http://www.mathworks.ch/ch/help/images/ref/fspecial.html */
double gaussian(double x, double y, double sigma){
	double res;
	res = exp(-(pow(x,2.0)+pow(y,2.0))/(2.0*pow(sigma,2.0)));
	return res;
}

/* normalize values */
void normalize(double **values, int size, double target){
	int x,y;
	double sum;
	sum=0.0;
	for(y=0; y<size; y++){
		for(x=0;x<size;x++){
			sum += values[y][x];
		}
	}
	for(y=0; y<size; y++){
		for(x=0; x<size; x++){
			values[y][x] /= sum;
			values[y][x] *= target;
		}
	}
}

int main(int argv, char** argc){
	int nW = 3;
	double dSigma = 1;
	double dMean = (nW-1.0)/2.0;
	double **kernel;
	double dSum;
	double dNormTarget = 1;
	int x, y;

	/* override default values */
	if(1 < argv){
		dSigma = atof(*(argc+1));
		if(2 < argv){
			nW = atoi(*(argc+2));
			dMean = ((double)nW-1.0)/2.0;
			if(3 < argv){
				dNormTarget = atof(*(argc+3));
			}
		}
	}

	/* title and debug informations */
	printf("Convolution matrix kernel\n");
	printf("argv: %i\n",argv);
	printf("Sigma: %g W: %i normalize to: %g\n",dSigma, nW, dNormTarget);

	/* dynamic memory allocation */
	kernel = malloc(nW * sizeof(double*));
	/* create array */
	if(NULL != kernel){
		for(x=0; x<nW; x++){
			(*(kernel+x)) = malloc(nW * sizeof(double));
			if(NULL == (*(kernel+x))){
				printf("malloc error: *kernel %i", x);
			}
		}
	} else {
		printf("malloc error: kernel\n");
	}

	/* calculate kernel and normalize to one */
	for(y=0; y<nW; ++y){
		for(x=0; x<nW; ++x){
			kernel[y][x] = gaussian((double)x-dMean,(double)y-dMean, dSigma);
		}
	}

	normalize(kernel, nW, dNormTarget);

	/* print kernel */
	dSum = 0;
	for(y=0; y<nW; y++){
		for(x=0; x<nW; x++){
			dSum += kernel[y][x];
			printf("%.4g,\t",kernel[y][x]);
		}
		printf("\n");
	}
	printf("sum: %g\n", dSum);

	/* cleanup memory */
	for(x=0; x<nW; x++){
		free(*(kernel+x));
	}
	free(kernel);
}
